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ABSTRACT: We discuss the Fermion sign problem and, by examining a very 
general Hubbard-Stratonovich (HS) transformation, argue that the sign problem 
cannot be solved with such methods. We propose a different kind of transforma- 
tion which, while not solving the sign problem, shows more detailed information 
about the system. With our transformation it is trivial to tell which auxiliary field 
configurations give a positive sign and which give a negative sign. We then discuss 
briefiy various properties of this transformation and construct a new algorithm 
which with one simulation gives results for a whole range of particle densities and 
Hubbard U values, positive and negative. Our approach is in excellent agreement 
with exact calculations. 
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The major obstacle facing numerical simulation of a large number of strongly interact- 
ing electrons is the so-called "fermion sign problem". ■'^'^'"^'^'^'^'^ This problem appears in 
many different guises and here we will focus on its form in the determinant algorithm for 
quantum Monte Carlo. Our method and conclusions are general but for clarity we will con- 
centrate on the Hubbard model, which is under active study as a model for metal-insulator 
transitions and high temperature superconductivity. We begin by reviewing the Hubbard- 
Stratonovich transformation and the resulting sign problem. The partition function of the 
Hubbard model in the grand canonical ensemble is given by 

Z = trie-^H), (1) 

H = -t J2 riWca(i) + 4(jK(z)) + C/J]K(z)-i)(n_(z)-^) 

<ij>,c7 i 

(2a) 

i 

= J2ct(i)kijC,{j) + uY,in+(i) - l){n-{{) - I). {2b) 
ij,a i 
The sum < ij > is over all pairs of nearest neighbor lattice sites, t is the hopping parameter, 

cj-(z) and Ccr{i) are creation and annihilation operators of electrons with spin a along the 

z axis at site i. n+ (n-) is the number operator for electrons with up (down) spins, 

P is the inverse temperature, and U is the coupling constant, which can be positive or 

negative. The matrix k in eq.(2b) contains the hopping term and the chemical potential. 

Through standard techniques one can transform the trace into a path integral (or sum) 

over the configurations of a c-number auxiliary field. ^ First we use the Trotter-Suzuki 

approximation^ to express Z as: Z = tr{e~'^^)^ ^ tr{e~'^^e~'^^)^ , where r = ^ <C 1, 

is the imaginary time step, and L is the number of such time steps. Now we use the 

HS transformation to decouple the quartic potential term, e~^^, into quadratics in the 



creation and annihilation operators: 
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-tU/4 

g-r[/(n+(i)-i)(n_(i)-i) = ^ L ^ ^-Xs{i,l){n^(i)-n_(i)) (3) 

s{i,l)=±l 



tU 

at each site i and time slice I. A is related to U via cosh{\) = e 2 . Here we wrote the 
discrete HS transformation,'^ but one can also write a continuous one. The trace over the 
fermion operators can now be taken since they only appear quadratically. This gives:^ 

Z= detM+detM-, (4) 

s{i,l)=±l 

where 

M^ = I + BIBI_^...B'[, (5) 

and 

I is a,V X V unit matrix, V is the spatial volume, v{l)ij = Sijs{i, I), where i runs from 1 
to V and / from 1 to L. The partition function is now written as a sum over c-numbers 
and can therefore be simulated on a computer. The sign problem arises because the 
determinants and their product, in eq(4), can be negative and thus cannot be used as the 
probability density in a Monte Carlo simulation. Instead, their absolute value is used, 
and the average of an observable. A, is then given by < ^ >=< Asgn >' / < sgn >', 
where <>' denotes averages with respect to the absolute value of the determinants. The 
average sign is thus defined as < sgn >'= Z/Z', where Z' is the partition function resulting 
from using the absolute value of the determinants. This approach works well for small U 
and relatively high temperatures. As the temperature decreases (/? increases), the average 
sign scales like'^'^ < sgn >~ e"*^^, where c is a constant. This makes low temperature 
simulations impractical because averages are obtained as the ratios of two very small 
numbers, each with a very large variance. When U < 0, a, similar procedure yields = 
M~. Consequently, the product of the two determinants in eq (4) becomes a square, 
and therefore positive semidefinite even though the determinant itself is still not positive 
semidefinite. Thus, there is no sign problem for the negative U Hubbard model (or other 
attractive interactions) . 
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The HS transformation discussed above is the most commonly used one, but it is 
only one of an infinite number of possible transformations. For example, other decoupling 
schemes were discussed in Refs.(5, 6) in the hope of finding a transformation which will 
solve the sign problem or at least decrease its severity. We will argue here that there exists 
no HS transformation that will eliminate the sign problem. We start by noting that the 
purpose of any HS transformation is to decouple the quartic fermionic interaction into 
quadratic terms which are coupled to the HS (auxiliary) field. This allows us to perform 
the trace in the partition function, giving two determinants. But, except for minor details, 
all transformations examined so far have yielded either a product of different determinants 
of the above form, or a single determinant^ which have always suffered severely from the 
sign problem for certain values of (3 and U . This suggests that one way to solve the sign 
problem in this approach is to obtain a square of a determinant. Is it, therefore, possible 
to generalize the above HS transformations such that the resulting partition function is a 
sum (or integral) over a square, and if not, why? In order to get (detM)^, the relative 
minus sign between n_ and n+ on the right hand side of eq (3) must become a plus, thus 
preserving the up-down symmetry of the Hamiltonian in the HS transformation. Let us 
therefore propose a general HS transformation 



where P{y) and y are real^ and arbitrary except for the constraints discussed below, and 
U is positive or negative. This transformation includes discrete transformations like eq 
(3), but is more general. If such a transformation were possible, the sign problem would 
be solved because of the resulting (detM)^. Since n± = 0, 1, we see that the conditions on 
P{y) and y are 






(8) 



< y >= e 2 , 



(9) 
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< >= 1, (10) 

where <> means an average with respect to the weight P{y). In general, the inequahty 
< y >^<< y^ > must be satisfied. Combining this with eqs (9,10) forces U to be negative: 
In other words, the inequahty cannot be satisfied for positive values of U. Therefore, there 
is no general HS transformation that is capable of giving the square of a determinant. 
Consequently, this approach to solving the sign problem fails. 

Implicit in the above argument is the positivity of P{y)- We will reserve the name 
"Hubbard-Stratonovich transformation" for these cases since all previous applications of 
the HS transformation assume such positivity. However, this argument is invalid if we 
allow P{y) to take negative values. This of course would not solve the sign problem, but it 
would give us a square of a determinant, with the result that the sign changes now come 
not from the determinants, but from P{y). One advantage of this is that contrary to the 
HS transformation where, in general, we do not know which auxiliary field configurations 
lead to minus signs because the structure of the determinants is too complicated, here the 
minus sign comes from P{y) which we know exactly. Therefore, we have complete prior 
knowledge of the sign of all the configurations. 

For example, a simple choice is 

P{y) = n(«'^^(^(^' - ^i) + - 2/2)), (11) 

i,l 

where i is the space index, and I is the time slice index, yi and y2 are the allowed values 
for the aiixiliary field y, a and b are parameters to be determined from the constraints 
eqs (8,9, 10). We chose two discrete values, yi and y2, for the auxiliary field, but we could 
have chosen any number of discrete variables, or any continuous distribution as long as we 
satisfy conditions eqs(8,9,10). Our motivation is simplicity. Applying the above constraints 
gives a + b= 1 and 

6=4^' (12) 
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g— = _L^1^_ (13) 
2/1 + 2/2 

Having chosen the form of P{y\ we still have the freedom of choosing the values of one of 
the two parameters 2/1, 2/2, the second being determined by eq(13). 

Now that we have decoupled the quartic fermion interaction into two quadratic terms 
with the same sign, the trace over the fermi operators can be performed as in the HS case 
giving for the partition function 

Z= a''^b''^{detMf (14) 

y{'i,i)=yi,y2 

= a^j:<i-T<(detMf>n2- (15) 
m 

M has the same form as M~, eq(5), and 

Bi = v{l)e-^K (16) 

with v{l)ij — 6ijy{i,l). ni (77.2) is the number of auxiliary spins with the value yi (1/2)) 
and N = ni + n2 is the total number of sites on the (d+l) dimensional lattice. <>n2 is 
an average over all configurations which have n2 spins equal to y2, whose number is the 
binomial coefficient C^. It is easy to show that for U < 0, both a and b are positive and 
therefore there is no sign problem, just like the usual HS transformations. When U > 0, a 
and b have opposite signs (we take b < 0) and thus the sign problem reappears. However, 
although the value of (detM)"^ (always positive) depends on both the relative number of 
spins with values yi and y2, and their configuration on the lattice, the prefactor aP'^b'"'^ 
depends only on the relative numbers. In particular, since the source of the sign problem 
in our formulation is the opposite sign of a and b, we have complete knowledge of all the 
configurations that change the sign: only configurations with an odd number of y2 spins 
lead to an odd exponent for b and thus a negative contribution. This complete character- 
ization of the negative configurations is to be contrasted with all the HS transformations 
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previously used where one knows very little about the kind of configurations that lead to 
minus signs. This vividly demonstrates, yet again, ^ that the sign problem in the deter- 
minant algorithm is not related to configurations where the electron "paths" exchange. 
It is merely an artifact of the transformation used to decouple the quartic terms in the 
Hamiltonian. Another property of our transformation is that it preserves the rotational 
spin symmetry of the original Hamiltonian. We can therefore take measurements along 
any spin direction, or along all three, thus reducing fluctuations. 

Not surprisingly, the behaviour of the average sign and the summand in eqs (14,15) 
depend on the values we choose for the auxiliary fields. Three choices are of particular 
interest to us. The first is to choose yi ~ 2/2- This choice is remarkable because the entire 
phase space is explored with a variable that fiuctuates very little. Consequently the values 
of the determinant and other observables fiuctuate little, and change very smoothly as 
more spins are fiipped. The disadvantage of this choice is that a/b ^ —1 (always keeping 
a + b = 1) as yi -^1/2, which means that the contributions of the negative configurations 
are of the same size as the positive ones. This makes numerical simulations hard, but may 
offer the possibility of studying the system semianalytically. 

Our second choice is to take 1/2 = 0. What is intriguing about this choice is that when 
enough auxiliary spins have the value y2 — 0, there will be many realizations where the 
auxiliary field for at least one entire time slice is zero, preventing the percolation of yi 
spins and resulting in M = / and detM = 1. When this happens the observables will also 
have a trivial value. So, when such configurations become important, there will be many 
instances where the observables are trivial and it would be interesting to study the effect 
such configurations might have on phase transitions. 

The third choice is to take y = ±1. The limit of eqs (12,13) as yi — > +1 and j/2 ~^ ~1 
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is 

a=l(l + e'f), (17) 

b=l{l-e'^). (18) 

Notice that with this choice, M is no longer a function of the coupling constant U, since the 
coupling constant can appear only in y, a and b, and here we fixed y = ±1. Furthermore, 
recall that, in the matrix M, fj, appears in the form e^l^I for each matrix, where I is the 
V xV identity. Consequently, the matrix M' = e~^'^(M — /) is a function of temperature, 
but not U or /i. This major simplification allows us to collect a large number of realizations 
of M', and then perform data analysis for any U and //. The /i dependence is obtained by 
multiplying M' by e^f^. The U dependence appears through the ratio b/a, eqs.(15,17,18), 
in a way similar to perturbation theory for r small. We can thus obtain results for a wide 
range of values of U (positive and negative) and from only one computer run. Currently, 
we generate the realizations of M' as follows: start with all y = +1 and successively flip 
randomly selected spins to ?/ = — 1, calculating M' after each spin flip. This is done up 
to some maximum number of flipped spins, much smaller than the total number of sites. 
The reason is that the expansion of Z in 77-2 (eq.l5), which is similar to a perturbation 
expansion for r small, converges quickly and can be truncated. Then we start all over, 
and in this way we generate an ensemble of realizations of flipped spins over which we can 
average. The addition of importance sampling could greatly increase the eflficiency of the 
algorithm. 

We tested this algorithm on a 2 x 2 lattice and compared with exact results. In flg.(l) 
we show a plot of < > versus U, with the crosses showing exact results. Fig. (2) shows a 
similar flgure for the ferromagnetic correlation function, 5(0, 0). We used 256 time slices in 
order to eliminate the finite time step errors for comparison with exact diagonalization. In 
general, the finite r errors are 0{t'^U), as in the usual determinant algorithm. Note that 
as we move away from half filling the errors increase appreciably for positive U because 
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the sign problem is appearing in force. We also see that for U <0, where there is no sign 
problem, errors are very small. Recall that all of the shown curves were obtained from the 
data of only one run, and that those data contain all the information needed to measure 
all the equal time correlation functions for positive and negative U and a wide range of 
II. One way to get better results for t/ > away from half filling is to increase statistics. 
Another more efficient way is to include importance sampling, and to correlate the update 
of positive and negative configurations since we know their signs. This would not change 
the average sign but would greatly decrease its variance and that of all measurements ^. 
This work is in progress. 
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Figure Captions 

Figure 1. The average density, < n-|- >, vs. U for fx = —0.1, —0.2, —0.3, —0.5 and — 1 
as labeled in the figure. The crosses show exact results. 

Figure 2. The ferromagnetic correlation function, S{0, 0), vs. U for the same values of 
the chemical potential as in figure 1. The crosses show exact results. 
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